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ABSTRACT 

Many global circulation models predict supersonic zonal winds and large ver- 
tical shears in the atmospheres of short-period jovian exoplanets. Using linear 
analysis and nonlinear local simulations, we investigate hydrodynamic dissipation 
. • mechanisms to balance the thermal acceleration of these winds. The adiabatic 

CLc Richardson criterion remains a good guide to linear stability, although thermal 

Q . diffusion allows some modes to violate it at very long wavelengths and very low 

growth rates. Nonlinearly, wind speeds saturate at Mach numbers ~ 2 and 
^. Richardson numbers ^1/4 for a broad range of plausible diffusivities and forcing 

strengths. Turbulence and vertical mixing, though accompanied by weak shocks, 
^ . dominate the dissipation, which appears to be the outcome of a recurrent Kelvin- 

ON , Helmholtz instability. An explicit shear viscosity, as well as thermal diffusivity, 

x/^ . is added to ZEUS to capture dissipation outside of shocks. The wind speed is 

^-^ not monotonic nor single valued for shear viscosities larger than about 10"'^ of 

the sound speed times the pressure scale height. Coarsening the numerical reso- 
O . lution can also increase the speed. Hence global simulations that are incapable 

of representing vertical turbulence and shocks, either because of reduced physics 
or because of limited resolution, may overestimate wind speeds. We recommend 
that such simulations include artificial dissipation terms to control the Mach and 
Richardson numbers and to capture mechanical dissipation as heat. 

Subject headings: binaries: close — hydrodynamics — instabilities — planetary systems- 
shock waves — turbulence 



1. Introduction 

Strongly irradiated extrasolar giant planets — "hot Jupiters" — are expected to rotate 
synchronously with their orbits but to have strong longitudinal winds that redistribute heat 
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from the day to the night side. The efficiency of redistribution is important for direct observ- 
ables including infrared phase curves, the depth of secondary echpses in transiting systems 
(i.e., echpses of the planet by its star), and planetary spectra. Tu rbulence associated with 



the winds may contribute to chemical mixing of the atmospher e (ISpiegel et al.l 12009), an d 



might even inject heat into the convective interior of the planet (iGuillot fc Showmanll2002l ). 
thereby perhaps explaining why some of these planets have lower densities than expected by 
standard evolutionary models at their estimated ages. 



In a previous paper (jGoodmanll2009l . hereafter Paper I), one of us has argued that these 
thermally driven winds can be understood as natural heat engines, which convert a fraction 
of the thermal power into mechanical work: namely, the work expended to accelerate the 
wind. As with any other heat engine, the continual production of work must be balanced 
by mechanical dissipation, else the kinetic energy in the winds would grow without bound. 
Paper I offered a brief discussion of possible dissipation mechanisms. Here we explore those 
mechanisms in greater detail. The goal is to lay the foundation for subgrid models of the 
dissipative process suitable for use by global circulation codes. 

The outline of this paper is as follows. §2]gives an overview of the dissipation mechanisms 
we consider, ^presents a linear analysis of hydrodjTiamic Kelvin-Helmholtz instabilities of 
thermally stratified shear flows. We show that such flows can be destabilized by thermal 
diffusivity even at Richardson numbers greater than the well-known critical value Ricrit = 
1/4. Under conditions relevant to exoplanet winds, however, we estimate that the associated 
growth rates are too slow and the longitudinal wavelengths too long to be important for 
dissipation, unless Ri < 1/4, which requires transsonic flow unless the vertical width of the 
shear layer is small compared to a pressure scale height. Hence shocks may be as important 
as shear-driven turbulence for dissipation. To investigate this, as described in ^ we have 
used the ZEUS code in two dimensions to simulate a part of the atmosphere with horizontal 
and vertical dimensions comparable to the pressure scale height. Thermal diffusion and 
viscosity have been added to the code, and thermal driving of the wind is simulated by 
a horizontal body force with nonzero curl. We study the velocity and dissipation rate of 
the wind in a statistical steady state as a function of the strength of the driving compared 
to the acceleration of gravity. §S] discusses our results in the context of previous work on 
winds in jovian planets both within and beyond the solar system. §H] summarizes our main 
conclusions. 
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2. Overview of dissipation mechanisms 

We begin by reviewing the physical conditions in those parts of the atmosphere where 
strong winds may occur. Because of the strong stellar irradiation, the temperature should 
be approximately constant with depth, 

where L^, is the stellar luminosity, a the orbital semi-major axis, and / a dimensionless 
factor of order unity that depends upon the local zenith angle of the star, the albedo, and 
the efficiency of lateral heat redistribution (/ = 1 for complete redistribution and negligible 
albedo). For a solar-mass primary, the orbital period is Porb = 3.66(a/lOi?0)^/^ d. The sound 
speed in an atmosphere dominated by molecular hydrogen is Cg ~ 2.75(T/1300K)^/^kms~^. 
A characteristic circulation timescale is 

where Rp is the planetary radius. Assuming that the planet rotates synchronously, by this 
definition tcirc is comparable to the rotation period, as it is on Earth. Coriolis effects are 
likely to be important but not as dominant as on a more rapidly rotating and colder planet 
such as Jupiter itself 

Several considerations point to transsonic wind speeds, as discussed in Paper I. First, 
absent dissipation, flow along isobars from the day to the night side, starting from rest, 
would accelerate to 

^ night 

where Az ~ Hp = c^/'yg is the range of depths over which the thermal forcing changes sign 
[see eq. (TT6|) ]. Indeed, repeated cycling between day and night could accelerate the wind 
indefinitely, were there no mechanical dissipation. Second, the strong entropy stratification 
tends to stabilize the flow against turbulent dissipation at small Mach numbers. And third, 
many numerical models of exoplanetary circulation do predict supersonic winds (see Paper 
I for references). 

Another direct consequence of the radiatively induced temperature ([1]) is that the ver- 
tical density proffie is approximately exponential, with a scale height Hp = k-oT / gmH2 that 
is small compared to the planetary radius: 

^P ^on.. in-3 /^ T \(Rp\ V^ 



^ ^ 2.9 X 10-^ -^ i:^ iL (3) 

Rp V1300K; \Rj) V^J, 
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where g = GMp/R^ is the surface gravity and gj = GMj/R'j^^^ ^ 2500 cms~^. 
The pressure at the infrared photosphere is 

Pph = 1^ - 1.7 (^^Y (^) mbar, (4) 

3Ka Vcm^g V \9jJ 

where Kr is the Rosseland mean opacity appropriate to the temperature ([I]) . A larger charac- 
teristic pressure is that at which the thermal time is comparable to tcirc- Taking tth = H^/x, 
where x = IGcT^/i^p^cp is the thermal diffusivity and cp ~ 3.5kB/mH2 the specific heat 
capacity, we estimate the latter pressure at 

.-™(^)"(i3^)"(i;) (!;)"■"*- '^' 

The Re ynolds number of thes e winds is enormous. The dynamic viscosity of molecular 
hydrogen is Jstiel fc Thodoslll963h pu = 1.9 x 10~^(T/1000K)°*'^gcm-i s-\ and therefore 



IT , V / T^ \ -U.iO / 



i?e = r^rz ~ 5 X lOM ^— ^— i^ . (6) 

V Vmbar/ VlSOOKy \gj J ^^ 

At pressures above a microbar at most, therefore, it is appropriate to regard the winds as 
inviscid. 

By other dimensionless measures, however, these flows are less ideal. At these temper- 
atures and pressures, thermal ionization of alkali metals yi elds magnetic Reynol d s num bers 



Rm = CsHp / ri ^ (1), where 7] is the magnetic diffusivity. iBatygin fc StevensonI (120 lOf ) and 



Perna et al.l ( l2010l ) have suggested that non-ideal MHD effects could exert a significant drag 
on the winds and perhaps even contribute to heating the convective interior in the presence 
of a background planetary magnetic field > 10 G. 

A dimensionless inverse measure of thermal (rather than magnetic) diffusion is the Peclet 
number 

Pe ^ ^ « 2 f ^^V (^^V^^VY^V (7) 

X ■ Vmbar/ Vcm^g-V V1300Ky \gjj ' ^ ' 

Thus, at the pressure (|5]) where thermal and advection timescales may be comparable, Pe ~ 
10^, with no direct dependence on Kr. 

Another important dimensionless quantity is the Richardson number, which character- 
izes the relative strength of stratification and shear. The local Richardson number at altitude 
z is Ri{z) = N"^ / [dU / dzY , where N{z) is the Brunt- Vaisala frequency [eq. fITT]) ]. In the limit 
of infinite Re and Pe, a necessary condition for incompressible shear-driven instability is that 
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Ri < 1/4 somewhere in the flow (JDrazin fc Reidlll98ll ). For a flow that changes smoothly 



by AU over a vertical distance Az, a characteristic value for Ri in a constant-temperature 
background is, using eq. (ITT]1 and noting Hp = cl/'jg, 

(7-1) /A^"' 



where M = 2AU/cs is the Mach number associated with the amplitude of the shear profile. 
Thus since 7 ~ 7/5 for an atmosphere dominated by molecular hydrogen, M > 0.9{Az/Hp) 
is required for instability. This linear stability criterion can be related to an energy principle 
and therefore should govern nonhnear instabilities also. However, when the stratification 
is primarily thermal rather than compositional, as we expect for the wind zone of these 
atmospheres, its stabilizing influence can be undercut by diffusion of heat. Thus instability 
may be possible at Ri > 1/4 if the Peclet number ([7]) is not too large. One of the goals of 
the present study is to quantify this statement for hot- Jupiter winds. 

From this survey of characteristic physical conditions and dimensionless numbers, three 
possible sources of mechanical dissipation for the winds come to the fore: (1) MHD effects, 
probably requiring B >10G] (2) shear-driven turbulence, requiring Ri < 1/4 and/or moder- 
ate Pe; and (3) shocks, requiring M > 1 as we have defined it. Anoth er possibility, raised in 



Pape r I, is the dou ble-diffusive Goldreich-Schubert-Fricke instability (JGoldreich fc Schubert 



19671 : iFrickd 119681 ). which we do not consider here for two reasons: first, although the hn- 
ear GSF instability is axisymmetric, we beheve that its saturation can be studied reliably 
only with high-resolution 3D simulations, whereas we hmit ourselves to 2D; and second, 
because the instability shuts off when (in the usual cylindrical coordinates) dv^/dz = and 
d{wv^y/dw > 0, which conditions are not inconsistent with strong vertical shears at equato- 
rial latitudes. Nonlinear MHD effects present even greater numerical challenges. Therefore, 
we concentrate on Richardson turbulence and shocks. 



3. Linear Analysis 

We are interested in the hydrodynamic stability of a horizontal shear flow. We first 
consider the Boussinesq limit in which we neglect density variations except where coupled 
to gravity. The equations of hydrodynamics become 

-— = -Vp + Oez + uV'^v, (9a) 

dt 

^ = -N^e, . V + xV'^, (9b) 

V-v = 0. (9c) 
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Here kinematic viscosity v is taken to be a constant, 



e^-g 



dlnT 



6T 

rji J rri 



(ideal gas) 



(10) 



is proportional to the fractional temperature variation, 



N' 



p 



dlnT\ dlnr/d2 
dlnp) g' dlnP/d2 



—)■ 



1-^9 



7 



Hr, 



(ideal gas & Tq =constant) (11) 



is the square of the Brunt- Vaisala frequency of the background hydrostatic state, and the 
pressure scale height 



H„ 



dlnp 



d^ 



^B^o 



^!^}^ = ^_ (12) 

^^'n^p9 19 

We assume the optically thick diffusion limit for radiative transfer in which x = IQaT'^/SnCpp'^. 
We further pick the opacity k, oc T^p'"^ so that x = constant. Note that in our energy equa- 
tion we assume no background temperature gradient. Other authors have studied the linear 
stabi l ity of the Boussinesq equations (see e.g..lGagdll973l : lDudislll974l : iMaslowe fc Thompson 
197ll : lKoppell[l963 : iMiller fc Gag3ll972l : lGagelll972l ). but they do not consider in detail the 
destabilizing effects of thermal diffusivity in the viscous, optically thick limit. We are inter- 
ested in the growth rates in the weak viscosity regime and their relevance to our numerical 
simulations. 



Assuming a plane-parallel background state with horizontal flow 

6*0 = = po, vo = U{z)e^, 



(13) 



we can expand all perturbations with the dependence exp{—iujt + ikxX + iky-y) to obtain the 
linearized 6*^ order system of equations 



'Iz 



ik'^a 

V 

io 

X 



k^- 



ikM'' 



V 



la 



vu+(2k'-—jv';,+ 



{k^ 



■Oi, 



k']e. 



X 



-vu 



(14a) 
(14b) 



where a 



u 



kxU, k = {kl + k'^Y^'^, and primes denote d/d 



z. 



We use a relaxation method to solve this system of equations, and we test our results 
against limiting cases of background velocity profile U and parameters z/, x, and A^^ for 
which there are simple analytic solutions. See the Appendix for the limits and analytical 
solutions. We can nondimensionalize our parameters by defining the Richardon number 
as Ri = N'^Hp/U^^^, Reynolds number Re = HpU/u, Peclet number Pe = HpU/x, and 
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Prandtl number Pr = u/x = Pe/Re. (In §2, we scaled Re and Pe by the sound speed 
rather than the flow velocity in the expectation of a Mach number of order unity, but that 
is not appropriate in the present Boussinesq context.) We give all our linear results in terms 
of these dimensionless numbers. 

For the problem of interest we specify boundary conditions for the system (fT^ as follows 
We impose Viz = aX z = ±Zmax so that the fluid perturbations do not penetrate the 
boundary. We further impose v"^ = —ikxv[^ — ikyv[ = 0, corresponding to vanishing 
viscous stress, and 6i = 0, corresponding to fixed temperature, at both boundaries. 

Consider now the limit z/ — )■ 0, x — ^ of eqns flT^ 

From WKB arguments, the time taken for a wavelike disturbance to reach z = from z = 
2;i 7^ is infinite if Ri > 1/4. Further, from energy arguments the potential energy barrier 
separating two fluid elemen ts at different heights exceeds the available kinetic energy that 



can be liberated if Ri > 1/4 (jChandrasekharll 1 96 ll ) . From these considerations, we expect the 



stratified shear flow to be stable for i?i > 1/4 in the inviscid case with no thermal diffusion. 
As we turn up the thermal diffusivity the potential energy barrier between fluid elements at 
different heights can be more easily overcome, and we expect instabilities to survive to higher 
Richardson number. We expect there to be a critical Richardson number, above which the 
stratified flow is stable and below which there exist unstable Kelvin-Helmholtz like modes. 
We refer to this critical solution as the marginally stable mode. 

For our stability analysis we pick background velocity profile U = UQta.nh.{z/Hp). 
Squire's theorem tells us that in the limit of vanishing thermal diffusivity the largest growth 



rates are obtained for kyHp = (JDrazin Sz Reidlll98ll ). so we might expect this result to gen- 



eralize to the case with thermal diffusion. In Fig [T] we show the marginally stable Richardson 
number as a function of k^Hp for kyHp = 0, .1, .2 and Pe = 10, Re = 10'^. As \kyHp\ increases, 
the critical Richardson number decreases for all k^Hp. Growth rates at fixed Richardson 
number correspondingly decrease. We set kyHp = hereafter because we are interested in 
the fastest growing modes. Note that we find two distinct instabilities at long and short 
wavelengths. The long wavelength instabilities have not converged with z^ax- The solid 
lines indicate solutions computed with z^ax = ^Hp, and the dashed line indicates a solution 
computed at Zmax = 5Hp. The long wavelength instability is driven to l onger wavelengths 



as Zmax increases. We run into similar issues encountered by iDudid (119731 ) in the long wave- 
length limit. We must integrate to larger and larger Zmax or use an appropriate asymptotic 
boundary condition as k^Hp — )■ to obtain better convergence. 

Although our code has not converged for long wavelengths, our results suggest that there 
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Fig. 1. — Marginally stable Richardson number as a function of k^Hp for Pe = 10, Re = 10^, 
and at a sequence of \kyHp\ = 0, .1, .2. As \kyHp\ increases the critical Richardson number 
decreases at all kxHp. The dashed curve shows a separate run at Zmax = 5Hp, indicating the 
long wavelength results have not converged. 
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may be a real inviscid instability in the limit k^Hp — )■ 0. This is demonstrated analytically 
in Appendix B. The dimensionless form of the equations of motion given in the Appendix 
suggest that the long-wavelength modes should be sensitive to the combination of parameters 
PeRi k z^^/Hp, and the analysis guarantees the existence of the long-wavelength modes only 
where this combination is <C 1; some version of the long- wavelength modes may exist more 
generally, but we have not proven it. So it is not surprising that our numerical results for 
these modes have not converged with respect to -Zmax- 

In Fig [2] we show the critical Richardson number as a function of k^Hp at Re = 10^ 
and at a sequence of Peclet numbers. As the Peclet number decreases below ~ 40, the long 
wavelength instability begins to strengthen and survive above Ri = 1/4. The dashed line is 
a run at higher z^ax = ^Hp and shows the convergence of short wavelength results. As the 
Peclet number decreases below ~ 2.5 the flow at short wavelengths begins to destabilize above 
Ri = 1/4. For high Reynolds number and high Peclet number our short wavelength results 
are in rough agreement with the inviscid results of Dudis (1974). There are quantitative 
differences due to our assumption of an isothermal background temperature profile. 
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Fig. 2. — Marginally stable Richardson number vs. wavenumber k^Hp, for fixed Reynolds 
number Re = 10^ and a sequence of Peclet numbers Pe = {v /x)R^- At Pe < 40, a long- 
wavelength instability that survives at Ri > 1/4 develops. At Pe < 2.5, the short-wavelength 
modes destabilize above Ri = 1/4. The dashed curve indicates a run at higher Zmax = 5Hp. 
The short-wavelength results have converged reasonably well, and for Pe ~ 10^ — 10^, the 
critical Richardson number remains fixed at Ri = 1/4. 
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in Fig|3]we show the growth rates of the modes with Ri = ^Ricnt as a function of kxH. 



p' 



for fixed Re = 10 and at a sequence of Peclet numbers. The long wavelength growth rates 
at k^Hp = 3.2/Re are well over two orders of magnitude smaller than shorter wavelength 
growth rates. This result agrees with our estimation that the inviscid instability in the 
limit kxHp — >■ has weak growth rate. For our numerical simulations we will be primarily 
interested in the more unstable short wavelength instabilities that have maximum growth 
at kxHp ~ 0.5. Note that -Ricrit increases with decreasing Pe, so we are computing growth 
rates at higher Richardson number for the lower Peclet number curves. At fixed Ri, growth 
rates increase with decreasing Pe, as we expect from simple heat diffusion arguments. 

We conclude from Figs [2] and [3] that the linear growth rates of the fastest growing mode 
are unchanged in the parameter regime we are interested in, Pe ~ 10^ — 10^. 
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Fig. 3. — Growth rates of modes with Ri = ^Ricrit as a function of kxHp, for fixed Re = 
10^ and Pe = 1.25,2.5,12.5,167. The long wavelength instability that survives to high 
Richardson number has growth rates over two orders of magnitude weaker than shorter 
wavelength growth rates. 
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Numerical Simulations 



We use ZEUS-2D v2.0 fIStone fc Normanlll992l ) to run hydrodynamical simulations of a 
horizontal wind profile in a local equatorial section of the atmosphere. Cartesian coordinates 
X and z correspond to longitude and altitude, respectively, each extending over a few pressure 
scale heights: typically 20ifp x bHp. Longitude x is periodic, though with a period much less 
than the true circumference of the planet. The simulations include vertical gravity, thermal 
diffusion of heat, shear viscosity, and a horizontal body force per unit volume pf{z)e.x [eqs. 
( IT7|) & ( TT9|) ] rather than thermal forcing. A few words of justification are in order. 

The pressure acceleration is nonpotential due to longitudinal entropy gradients; the 
component of its curl normal to the equatorial plane is, in spherical coordinates r6(j), 



Be ■ (Vp-^ X Vp) 



1 fdp\ 


'dTdp 


dpdT 


_d(j) dr 


d(j) dr 



g I'dlnT 

r 



(16) 



in which the final form uses the ideal-gas law and assumes approximate vertical hydrostatic 
equilibrium. This curl must vanish in the (nearly isentropic) convection zone. Furthermore, 
since thermal forcing by itself does not add any net momentum to the fiowj^ the nonpotential 
longitudinal force must reverse sign with depth within the radiative atmosphere. In fact one 
expects this nonpotential "battery" to be concentrated at pressures < pth [dS])], which is quite 
shallow compared to the radiative-convective interface, typically ~kbar. The nonpotential 
force must also reverse sign at least twice in longitude because (fT6l) is periodic in 0, but the 
horizontal scale of variation of this force is probably small compared to its vertical scale since 
the radiative atmosphere is thin. Hence the profile ( IT9l) . which is horizontally constant and 
extends over ~ 2Hp vertically, reversing sign within this depth, is a plausible representation 
of the horizontal forcing within our local computational domain. 



4.1. Setup 



The source terms for momentum and energy in ZEUS are modified by the addition of 
the terms within square brackets below: 



dv 



-Vp - pge, - V • Q + [pfiz)e^ + V • cr'] 



(17) 



^ Arras fc Socratea ( 20101 ) have suggested, however, that in conjunction with the tidal field of the host 
star, thermal forcing may lead to net torques on the planet. 
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and 

^ ) = _pv. t; - Q • Vi; + [V ■ {xpCpVT) + ct' • Vv\ , {11 



in which <t' is the viscous shear tensor, with components 

dvi dvk 2 dvi 



where v and g are taken to be constants. ZEUS' standard von Neumann & Richtmyer 
artificial viscosity tensor Q is apphed with shocks spread over ^ 2 zones. The thermal 
diffusivity x is again given in the optically thick limit by x = 16o"T^/3/tCpP^ = const; 
constant x would correspond therefore to k oc T^p~'^. The horizontal forcing profile is 

f{z) = eg [2sech.^ (z/Hp) tanh{z/Hp) + asecli^iz/Hp)] , (19) 

where Hp is the pressure scale height, a is a parameter that we adjust to ensure zero net 
momentum input on our grid, and e <^ 1 determines the overal strength of the forcing. This 
form is motivated by the viscous acceleration on a fiuid with hydrostatic density profile and 
horizontal velocity profile given by u^. = UQta.nh.{z/Hp). The viscous acceleration in this case 
is given by eg = vlIo/Hp and a = 1. We use an FTCS differencing scheme for our modified 
source terms, and the ZEUS timestep constraint must satisfy 5t < Sty = [AxY/Au and 
6t < 6t^ = (Ax)^/4x for stability. We are interested in Pr < 1, so the latter constraint is a 
tighter restriction on 6t. But St^/6tcFL ~ P^/4:Np, where Np is the number of grid points 
per pressure scale height and Pe is the Peclet number in the fiow. Dimensionless parameters 
for our numerical results are scaled to the sound speed Cg as in ^ Our simulations generally 
satisfy Pe ~ 10 - 100 and Np ~ 0(10), so 5tJ5tcFL ~ 0(1) and the diffusion and CFL 
timesteps are usually comparable. 

As stated, x is periodic. At z = z-min and Zmax, we impose reflecting conditions for 
density and velocity and constant temperature T = T^aii- We initialize the density to 
p(x, z) = po and the internal energy density to e{x, z) = cq, allowing the flow to settle 
under the influence of gravity. The velocity is initialized to t7 = Uotanh.{z/Hp)ex. The wall 
temperature is set to the initial grid temperature, T^aii = ^o = (eo/Po)/^"^p(7 ~ ^)/kB- We 
pick po = I5 Co = .01, 7 = 1.4, and g = .004. Our simulations are run over 5 vertical pressure 
scale heights, from z = —2.5Hp to z = 2.5Hp. The density at the base of our simulation if 
it were in hydrostatic equilibrium is then pbase = 5po- The code units for our problem are 
given in table [U We nondimensionalize parameters according to these units. 

We run our simulations over 20 scale heights in the horizontal direction. We are phys- 
ically interested in Pe~^ ^ 10^^ — 10~^, Re~^ — )■ 0, but the viscosity must be large enough 
to maintain energy conservation, and e < 10"^. In the absence of a shear viscosity the grid 
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loses energy to numerical dissipation at the grid scale. We must pick v large enough such 
that the viscous dissipation length l^ ~ (|(9jfi|/z/)~^/^ is larger than the grid scale, and the 
dissipation of solenoidal turbulence is converted to internal energy. Our standard runs have 
A'p = 10 grid points per pressure scale height, but for certain runs which conserve energy 
less well we use up to N^ = 30 grid points per pressure scale height. 



4.2. Tests 

We run two simple test problems to verify proper implementation of the horizontal 
forcing and viscosity. Using a simplified shear viscosity dv/dt = z/V^f and no viscous heating, 
there is an exact solution for the viscous decay of initial velocity profile given by U{z, 0) = 
sign(2;)e^. The velocity is given at later times by U{z,t) = eTi{z/2\/ui)e.x. Again using this 
simplified shear viscosity and horizontal forcing given by f{z) = e2sech.'^{z/Hp)taiih.{z/Hp), 
there is an exact laminar solution for the velocity profile. The laminar solution satisfies 
= f{z) + uU"{z), where U{z) = {eHp/h')taiih.{z/Hp). We conclude from these two tests 
that there is negligible numerical dissipation in the laminar regime. We can also use the latter 
to estimate the wind speed in the laminar regime when the viscous heating is implemented. 

We are interested in the behavior of the solutions to eqns [17] and [I8] at different forcing 
amplitudes. We check that the parameter a has a second order effect on the forcing and 
that we can roughly characterize the amplitude of the forcing by the single parameter e. 



4.3. Diagnostics of dissipation 

We must also maintain energy conservation in our simulations. We are continually 
injecting energy onto the grid via the horizontal forcing, and to reach a statistical steady 
state this energy must first be converted into internal energy and then allowed to diffuse 
through the boundary walls. Because ZEUS is not based on a conservative form of the energy 
equation, numerical diffusion dissipates kinetic energy on the grid scale without converting it 

Table 1. Code Units. 



Unit 


Length 




Time Density 


Acceleration 


Value 


Hp = l 


t = 


= Hp/Cs = 13.36 Phase = 5 


g = .004 
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to internal energy; to avoid this, it is necessary to add an explicit shear viscosity (in addition 
to the artificial bulk viscosity used to mediate shocks). 



We compute a fractional energy error 

-C/ -C/in ~ -C/kin ~ -C/int ~ -C/prav ~ -C/n 



where 



Eu 



E, 



out 



^grav -'^out 



Ew 



pvxf\z)dxdz 



-pCpX 



or 

dz 



E. 



E- 



kin 



dt 
d 



1 „„,2 
2 



dx -Egrav = ~f 1 1 P9^ dxdz 



(20) 



pv^dxdz E'int 



d 
dt 



edxdz 



We briefiy mention here a number of strategies that we use to improve energy conser- 
vation. The fractional energy error will tend to decrease if we keep the amplitude of the 
forcing function lower so that there is less energy injection onto the grid. In a similar vein, 
if we localize the forcing function to the shear layer as opposed to distributing it over the 
entire vertical extent of the grid then there will be less energy injection and better energy 
conservation. Other strategies include increasing the kinematic viscosity or making the grid 
resolution finer in order to capture kinetic energy dissipation above the grid scale via the 
explicit viscosity. Finally, we can increase the thermal diffusivity so as to allow excess in- 
ternal energy to diffuse off the grid. These considerations must be balanced against keeping 
the computation time reasonably low and keeping the solutions physically relevant to the 
problems of interest. 

We consider that a statistical steady state is reached when {Eint)t, (-E'kin)*, and (-Egrav)t 
are all <^ Ein, and when the Mach number [eq. ( I2^ ] and the rms vertical velocity {v1)x,y are 
constant for at least a horizontal sound crossing time, 20Hp/cs. In steady state, E/Ei^ ~ 
(iJin — -Eout)/-£'in- We maintain the fractional energy error to within ten percent. 

Dissipation of kinetic energy involves production of entropy. If s is entropy per unit 
mass, then 

Ds -Q : Vv + a' : Vv ,„, ^,2 „ , „, rr.^ 

p— = + xcpp I V In T| + V-(xCppV In T). 

When integrated over our computational domain, the final term is — -Eout/^waii, -E'out being 
the rate of escape of heat. So 

^out = T^aii / [T-\-Q :Vv + (7' : Vv) + xcpp| VlnTp] d^x + ^ ffpsdxdy. (21) 
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In statistical steady state, the time average of the last term vanishes. 

While conversion of mechanical energy to heat ultimately depends upon the "micro- 
scopic" processes in eq. fl^Tl) . we gain insight into the macroscopic processes that remove 
kinetic energy from the large-scale flow — shocks and turbulence — by monitoring the energy 
change in each operator split source step. The price paid for this approach is having to deal 
with a mixture of reversible and irreversible effects. Under adiabatic conditions, for exam- 
ple, the compressive heating term Epdv = f — P V- v (Px would be completely reversible, but 
shocks and thermal diffusion make {Epjy) 7^ even in steady state. There is also a change in 
gravitational potential energy -Egrav = / Vzpgd^x associated with the vertical gravity source 
step. In many of our runs we find that Epdv and -Egrav have large oscillations that nearly 
cancel. Using integration by parts, one can show that 



E,^ 



pdV 



+ 4rav = / / ^ ■ [P^iaz] + Vp] dxdz . (22) 



Thus the near cancellation suggests that our flows are close to vertical hydrostatic equilibrium 
despite large vertical motions. We find it easiest to monitor the combination Ep^v + -E'grav 
The time average of this term is largely the work done to mix the fluid vertically; ab- 
sent thermal diffusion, this could lead to an almost adiabatic temperature-pressure profile, 
dlnT/dz ^ —(7 — 1)(7/Cg, but in steady state diffusion tends to restore the stratification. 
Thus vertical mixing acts as a refrigerator — a heat engine in reverse — absorbing mechanical 
work to drive heat up the temperature gradient. 

Shock dissipation is dominated by the artificial bulk viscosity, -Eabv = / —Q • 'Vvd'^x 
but receives small contributions also from the shear stress and from the Epdv term. 

Since we find that shocks are weak, we can estimate the turbulent energy dissipation 
rate as 

{Eturh)t ~ {EpX^.v + Egri,^)t + {E,y)t — (-Elam)*, (23) 

where Ey = j cr' • Vv and -Eiam = / v p{dvx / dz^d'^x. Assuming an effective turbulent vis- 
cosity acting on the time- averaged mean shear profile, we estimate the turbulent dissipation 
rate as (-Eturb)t = (/ J-'tp{,dvx / dz)"^ d'^x) t. The effective shear viscosity due to turbulence in a 
flow that does not resolve the turbulent energy cascade is then u^ = i^{EtuTh)t/{Eia.m)t- 



4.4. Results 

Figs Hl[5| and [6] respectively give the run of Mach number, defined in eq. fl2^ . rms 
vertical velocity, and fractional energy error versus u/cgHp for x/csHp = 0.094, e = 0.1. The 
Prandtl number is less than one for all of these runs. Note the presence of four distinct 
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solutions at v/csHp > 10~^. There is a fast solution, denoted by filled triangles, a medium 
velocity solution, denoted by filled squares, and two slow solutions, denoted by filled circles 
and open stars. Where solution branches overlap, we trace each branch by starting from a 
representative member and then gradually varying parameters. We use a similar parameter 
deformation to verify that each branch is limited to a certain range of u/cgHp as indicated 
in Figs mH 
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Fig. 4. — Wind amplitude as a function of scaled viscosity u/cgHp for e = 0.1, x/cgHp = 
0.094. We find four distinct solutions in the range log^g ^/(^sHp = — 3 to — 1.5, but a unique 
solution, with Mach number Ai ~ 2.3, in the high-i?e limit. 

The solutions differ in the importance of laminar, turbulent, shock, and irreversible 
compressional dissipation. The fast solutions (triangles) have high laminar energy dissipa- 
tion. Fig |5] indicates that there is a significant vertical component to the velocity field, 
but the bulk of the viscous dissipation is laminar rather than turbulent. This is reflected 
by wavelike oscillations in the flow as opposed to turbulent Kelvin-Helmholtz rolls. The 
medium and slow B solutions have stronger irreversible compressional heating and shocks, 
but they still have relatively high laminar dissipation. For the planetary wind problem of 
interest to us this laminar dissipation is unphysical, and as we increase the Reynolds number 
these solutions disappear. We conclude that these unphysical solutions are supported by 
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Fig. 5. — Scaled vertical velocity fluctuations as a function of u/cgHp for e = 0.1, x/^sHp 
0.094. 
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Fig. 6. — Fractional energy error as a function of u/cgHp for e = 0.1, x/c^ifp = 0.094. 
The extra two lower filled circles at log^Q u/cgHp = —3.03 demonstrate the convergence of 
our slow A solutions. The energy error decreases with increased resolution, but the gross 
properties of the flow remain unchanged. 
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the shear viscosity. Numerical diffusion can also support these unphysical solutions. At a 
grid resolution of A'p = 3 grid points per pressure scale height and parameters from Fig 
m we find unphysical solutions even at v/cgHp = 0. We estimate that it is necessary to 
have Np > 0(10) in order to enter the high-i?e regime where there exists a single, physical 
solution. This unique solution is the slow A solution in Figs H] - [61 It is characterized by 
weak laminar dissipation and strong irreversible compressional and shock dissipation. We 
check the convergence of this solution with runs at 1.5 and 2 times the standard resolution. 
Extra filled circles at log^Q u/cgHp = —3.03 in figs H] - El indicate higher resolution runs. The 
wind speeds and rms vertical velocities overlap, but the energy errors progress downwards 
at higher resolution. 

We summarize these results in Table [2] by providing quantitative diagnostics for each so- 
lution type. We give the shear viscosity, turbulent viscosity, relative importance of artificial 
bulk viscosity, relative importance of turbulent dissipation, Mach number, and Richardson 
number for a representative solution of each type. The laminar and turbulent energy dissipa- 
tion rates are directly proportional to the quoted shear and turbulent viscosities, as described 
in §4.3[ At this stage we are only providing rough diagnostics to characterize each solution 
type, but we infer that shocks contribute weakly to the irreversible compression because the 
artificial-bulk-viscosity dissipations are relatively small. The Mach number is computed as 



M 



max 

2 



J p{x,z)v^{z)dx 
Cs,wi,iiJ p{x,z)dx 



— mm 



J p{x,z)v^{z)dx 
Cs,wi,iiJ p{x,z)dx 



(24) 



We compute the minimum Richardson number in our flows via 



Ri 



N^ 



mm 



(dU/dz) 



mm 



dlnp 
dz 



fdU 
\d~z 



-2 



^ V 



ad 



rflnT 
d\iip 



where overbar denotes horizontal average. The fastest and medium solutions are convective 
near the top wall. This is due to the higher wind speed and energy injection, which cause a 
steeper temperature gradient there. 

Fig. m shows that the fastest triangle solution in the laminar regime has wind speed 
scaling roughly as 1/z/, as we expect from our tests in §4.21 Further, the perfectly laminar 
solution at v/csHp = —1.03 has Ri = 0.36, and we check that it should be linearly stable as 
well. From our runs at lower forcing we determine that for the fastest (triangle) solution, the 
viscosity that maximizes the wind speed scales in proportion to the forcing: Vpeak/csHp oc e. 
The peak wind itself decreases only marginally. We can understand this behavior by assuming 
a steady-state balance between laminar dissipation of the mean flow and the forcing flTOj) . 
These solutions are able to reach such a high wind speed because they are nearly laminar, 
until turbulent dissipation kicks in at low v. 
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We map the solution space at lower forcing and thermal diffusivity and verify that there 
are four different solutions in the regime Pr < 1. Specifically, we looked at x/(^sHp = 
0.094 and e = .025, .05, and 0.1, and also at x/'^sHp = 0.047, e = 0.1. These parameters 
correspond to the four values in the first three rows of table |31 discussed in greater detail 
below. There are small quantitative differences in the results at x/^sHp = 0.047, e = 0.1. 
The fastest solution in particular has peak horizontally averaged Mach number ~ 15% higher. 
From our linear analysis we expect little difference in the results for any Peclet number in the 
range Pe > 10. In our simulations we are continually injecting energy onto our grid though, 
and the diffusivity affects how this energy is redistributed. We attribute the differences in 
the nonlinear behavior for Pe ~ 0(10) to the important dynamical effect of this internal 
energy redistribution. 

The main takeaway point from these runs at different forcing and thermal diffusivity 
is that in the high-Reynolds-number limit there is a unique solution, for which the direct 
viscous dissipation of the mean flow is negligible compared to turbulent and shock dissipation. 
Further, if the grid resolution is insufficiently low, with Np < 0(3), we may lock onto an 
artificially fast solution, even at u/cgHp = 0. We check the latter result for the range of 
parameters we studied and not just for those relevant to Fig. |H 

We now attempt to extract an effective turbulent viscosity from the physically relevant 
slow A solution. First note that the energy errors in Fig |6] are largest for the slow A solutions. 
In order to decrease the numerical dissipation we can either increase the shear viscosity 
v/cgHp or we can make the grid resolution finer. We cannot increase the shear viscosity 
indefinitely because we will enter the multiple-solution regime (Fig. |4]) or the laminar viscosity 
will become large enough to modify the properties of the solution. We must balance this 
against the increased computational cost of running the simulations at higher resolution. To 

Table 2. Viscosity, turbulent viscosity, relative importance of artificial bulk viscosity, 
relative importance of turbulent dissipation, Mach number, and Richardson number for a 

characteristic solution of each type in Fig HI 



Solution type 


log 


wu/csHp 


J^t/csHp 


{Eahv)t/{Ein)t 


{Etnih)t/{Ein)t 


M 


Ri 


Fastest (triangle) 




-2.17 


0.006 




0.07 




0.43 


7.6 


<0 


Medium (square) 




-2.57 


0.010 




0.22 




0.60 


5.5 


<0 


Slow A (circle) 




-3.03 


0.015 




0.13 




0.79 


2.3 


0.033 


Slow B (star) 




-2.40 


0.017 




0.12 




0.67 


2.9 


0.024 
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calculate the turbulent viscosity we use sufficiently high grid resolution and shear viscosity 
such that the numerical energy dissipation rate is < 5% of the energy injection rate, but 
shear viscosity is low enough such that the gross properties of the solution are not modified. 

In Tables [3] and H] we give the Mach number and turbulent viscosity for a series of slow 
A solutions at different horizontal forcing amplitudes and thermal diffusivities. The exact 
value of the explicit viscosity we use is unimportant, so long as the viscous dissipation length 
is larger than the grid scale and the viscosity is small enough to leave the gross properties 
of the slow A solution unmodified. In general the runs satisfy Re > 10^. The runs at 
x/csHp = 0.094 and e = 0.025, 0.01 are done at the standard resolution of Np = 10 grid 
points per pressure scale height. The rest are done at twice the standard resolution (i.e., 
twice as many points in each dimension), with the exception of the run at x/^sHp = 0.013, 
e = 0.01, which is done at three times the standard resolution, and the run at x/csHp = 0.094, 
e = 0.005, which is done at 1.5 times the standard resolution. We emphasize that these values 
of Np indicate the grid resolution required at each set of parameter values to conserve energy 
within 5%, not the resolution required to find the slow A solution at v/cgHp = 0. At the 
standard resolution of Np = 10, we still find the slow A solution at u/csHp = for all 
parameters surveyed, but the fractional energy error may be as large as 0.3 for the lowest 
diffusivity runs. 

The Mach number decreases gradually with forcing amplitude e for e = 0.1 down to 
e = 0.01, and is roughly independent of e for e = 0.005 — 0.01. It is also independent of 
the diffusivity strength for Pe ~ 10 — 100, at all forcing amplitudes surveyed. The latter is 
consistent with our linear analysis, and we attribute this result to the lower energy injection 
rate as compared to the faster solutions from Fig. HI There is less thermal redestribution of 
heat in the slow A solutions. 

Fig. [7] shows the Mach number as a function of time for the run at e = 0.01 and 
x/cgHp = 0.013. The run has reached a statistical steady state after time ~ lOHp/cg. The 
steady state is characterized by a drive-up phase in which the Mach number gradually rises, 
and a short slowdown phase in which the Mach number rapidly drops. In Fig. [8] we show the 
relative contributions of dissipation due to artificial bulk viscosity and turbulence, as defined 
in §4.31 Both dissipation terms are confined to the slowdown phases. Shocks appear briefly 
during slowdown, but the majority of the kinetic energy dissipation occurs in the absence 
of shocks. We thus conclude that shocks do not signiflcantly contribute to the irreversible 
compressional heating. 

We now characterize the nature of the sharp slowdown phases that limit the wind speed 
in our flow. The drive-up phase is relatively smooth with small changes in the state variables 
across the flow. Fig. [9] shows the velocity vector fleld with temperature contours from T„aii 
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Table 3. Mach number for a sequence of slow A solutions at different diffusivity and 

horizontal forcing. 



e 


.094 


x/csHp 
.047 


.013 


.1 


2.3 


2.3 




.05 


1.9 






.025 


1.7 






.01 


1.7 


1.6 


1.7 


.005 


1.5 


1.7 


1.7 





Table 4. Turbulent viscosity v^/csHp for a sequence of slow A solutions at different 

diffusivity and horizontal forcing. 



e 


.094 


x/csHp 
.047 


.013 


.1 


0.015 


0.017 




.05 


0.010 






.025 


0.0061 






.01 


0.0024 


0.0023 


0.0024 


.005 


0.0012 


0.0011 


0.0010 
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to Tmax = l-07Tmax foi the Starred point from Fig. [71 At a short time before reaching 
the peak wind speed, t ~ 0.2Hp/cs, a vigorous instabihty with 3 wavelengths in our box 
develops. Fig. [TU] shows the velocity vector field with temperature contours from T^aii to 
Traax = l-46Tinax for the boxed point from Fig. [71 The unstable mode is clearly visible 
in the temperature profile roughly halfway through the slowdown phase. Note the shocks 
at z/Hp < —1 and x/Hp = —4.5, 2, 8.5. Following our previous determination that the 
shocks are weak, we conclude that turbulent mixing associated with this unstable mode is 
the primary agent for dissipating kinetic energy and reducing the wind speed. 

The dimensions of our box enforce a periodicity under the transformation x — )■ x + 20Hp. 
If the unstable mode is an ordinary Kelvin-Helmholtz mode, then from our linear analysis 
the only modes that can grow have either 1,2, or 3 wavelengths in the box. Apparently the 
mode with 3 wavelengths and wavenumber k^Hp ^ .94 is the fastest growing of the three. 
The Richardson numbers satisfy Ri < 1/4. We verify that the growth of this instability 
and its effect on the wind speed are independent of the dimensions of our domain. We 
increase the horizontal extent of our box from L^ = 20Hp to L^ = 30Hp and L^ = AOHp 
and verify the presence of this instability. At L^ = 30Hp the fastest growing modes have 
3 or 4 wavelengths in the box, and at L^ = AOHp they have 4 or 5 wavelengths. These 
wave numbers are roughly consistent with our linear analysis. More importantly though, the 
time-averaged Mach number of the flow remains unchanged. 



5. Discussion 

Ultimately we would like to provide a dissipation prescription for global simulations 
that do not resolve the turbulent nature of the wind profile. We first note that from Fig 
[71 we find a horizontally and time averaged Mach number A^ ~ 1.7 and peak horizontally 
averaged Mach number A^ ~ 2.7. The peak occurs just before the rapid growth of instabil- 
ities, and so the latter corresponds to the absolute peak Mach number in the flow. Global 
simulations are often more strongly supersonic. For example, the model of HD 209458b by 



Cooper &: Showmanl (|2005| ) finds fmax ~ 9kms ^ at the 2.5 mbar level (see Paper I for a 



survey of the calculations prior to 2009). 

Most of these studies of global circulation employ reduced equations of motion based on 
the approximation of vertical hydrostatic equilibrium: shallow-water, equivalent barotropic, 
or primitive equations. This is done in order to cope efficiently with the large ratio of 
horizontal to vertical scales [eq. ([3])], but the approximations used also have the effect of 
filtering out sound waves and shocks (even horizontally propagating ones). Such codes cannot 
directly represent the shocks and turbulent vertical motions studied here, although they can 
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Fig. 7. — Horizontally averaged Macli number as a function of time for the slow A solution 
at x/csHp = 0.013, e = 0.01. The steady state after time ~ lOHp/cs alternates between 
phases of gradual acceleration, and sharp drops of Mach number due to a violent instability. 
The star and box denote points at which we illustrate the vector field in figs M and [TO] 
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Fig. 8. — Dissipation fractions via artificial bulk viscosity (abv) and turbulence at e = 0.01 
and the sequence of x/csHp from Table [31 The turbulent dissipation is largely irreversible 
compression, and it dominates over the shock dissipation. 
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Fig. 9. — Velocity field for the slow A solution at xl^^sHp = 0.013, e = 0.01 during the drive 
up phase (starred point in Fig. [9]). Equally spaced temperature contours are overlaid from 
T^aii to Tmax = l-OTT^aii- The minimum Richardson number is Ri = 0.11. 
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Fig. 10. — Like Fig. IH] but during the unstable phase (boxed point in FiglH]), and with 
contours from T„aii to T^ax = l-46Twaii. Minimum Ri = 0.04. Periodic shocks occur at 
z/Hp < — 1, but the bulk of the heating is due to irreversible compression. 
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represent horizontal shea rs (which we do not) and hydrauhc jumps (e.g lShowman et aLll2009l : 



Rauscher fc Menoull2009l ). which bear some resemblance to shocks but do not conserve energy 
and may not be triggered in the same circumstances. Therefore no particular reason exists 
as to why our results for the wind speeds should agree quantitatively with those obtained 
from the primitive equations, for example. 



Comp arison with the calcu l ations carried out by lBurkert et al.l ( l2005l ) , iDobbs-Dixon fc Lin 



( I2OO8I ) and iDobbs-Dixon et al.l (J20ld . herafter DCL) is potentially more revealing, because 
these calculations, like ours, use the full compressible hydrodynamic equations. The most 
recent of these studies included an explicit shear viscosity. It was found that the viscosity 
had little effect unless at least 10^°cm^s~^, above which the speed decreased monotonically. 
For reference, v = 10^°cm^s~^ corresponds to Re ~ 600 referred to the sound speed and 
scale height. We find, however, that the wind speed is not monotonic or even single-valued 
in the viscosity (Fig. (j4])). 

We would like to understand what accounts for the difference between our results and 
those of DCL. Compared to our body force and constant diffusivity, the thermal forcing and 
radiative transfer implemented by DCL is much more sophisticated and realistic, but it is 
hard to see what this has to do with the viscous effects, especially since the behavior seen in 
Fig. (jl]) exists for all choices of forcing and diffusivity that we have tried. Furthermore, our 
direct forcing is usually stronger than DCL's thermal forcing since, for example, a horizontal 
acceleration of / = 0.025(7 would result in a Mach number y^27rR^/cs > 7 after going 
half way around the planet without dissipation. And yet our Mach numbers are at most 
comparable to those of DCL, the latter peaking at ~ 3 for low viscosity. 

More likely the difference in results has to do either with the difference in the dimen- 
sionality and geometry of the calculations (3D and global for DCL; 2D and local for us), or 
with spatial resolution. It is well known that subsonic turbulence behaves differently in two 
and three dimensions; however, since two dimensions usually favor energy on large scales, 
our wind speeds might have been expected to exceed DCL's if dimensionality were the most 
important factor. Concerning numerical resolution, it appears that DCL's standard value 
is < 4 cells per scale height, compared to our 10. Since we and DCL use similar numerical 
methods, their calculations are probably more diffusive than ours. When we use a resolution 
comparable to theirs, our wind speeds actually increase for the same forcing ( §4.4p . The 
fact that we and they find comparable Mach numbers at our highest respective Reynolds 
numbers suggests that our stronger forcing is partly offset by our finer resolution. 

A limitation of our local, two-dimensional models is the somewhat arbitrary assumption 
that the nonpotential horizontal force is compensated over a depth range Az ~ 2Hp. Other 
things being equal, one would expect the flow speed to scale as (Az)^/^ at constant forcing 
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amplitude if Richardson number rather than Mach number is the controlhng parameter. 
Inspection of DCL's figures suggests that the horizontal speed scales approximately linearly in 
Inp over the range 10~^-10°bar. The sign change is not always obvious in their velocity plots, 
although they infer the existence of return flow from their temperature profiles. However, 
since the momentum density frvx oc Vx{z) exp{—z/Hp), this is much better localized around 
its sign change than is Vx itself. 

We can use our turbulent viscosities to estimate a turbulent diffusi vity for composition , 



Ki, as may be required to keep the upper atmosphere turbulently mixed (ISpiegel et al.ll2009l ). 



The Schmidt number, define d by Sct = Ut/Kt, in general has a complicated dependence 



on flow properties (see e.g., iHuq fc StewartI l2008l ). For Ri < 0(10 ^), however, Sct is 



independent of Ri and of order unity. The data in Table H] then gives us mass diffusivities 
pertaining to the sh ear layers in Hot Jupiter atmospheres: Kt > 10~^CsHp ~ 10^°cm^s~^. 



Spiegel et al.l (120091 ) estimated that Kt ~ 10^ — 10^^ cm^ s~^ is required, depending on grain 
size, to maintain TiO as an optical absorber at p < 1 mbar. Our values are toward the 
upper end of this range, but only in the shear layer (~ 2Hp). To estimate Kt elsewhere, we 
would need to estimate how much of the mechanical power dissipated in the shear layer can 
propagate upwards and downwards in the form of internal waves or shocks. This we have 
not attempted. 



6. Conclusions 

We have investigated the balance between horizontal forcing and mechanical dissipation 
in model shear flows designed to imitate thermally driven winds of jovian exoplanets. Our 
main conclusions are as follows. 

1. Following Paper I, some form of mechanical dissipation is required to offset the produc- 
tion of mechanical energy by the longitudinal entropy gradient. Plausible dissipation 
mechanisms include shocks, turbulence, and MHD drag. 

2. Entropy stratiflcation has a stabilizing influence on Kelvin- Helmholtz (KH) modes, 
but thermal diffusivity (x) tends to undercut this influence to a degree that can be 
quantifled by the Peclet number Pe = CsHp/x based on the sound speed (cg) and 
pressure scale height [Hp). We estimate Pe ~ O(IO^) at the altitude where the thermal 
timescale is comparable to the circumferential flow time at Mach 1, but Pe ^ 0(1) 
near the infrared photosphere. 

3. Linear analysis indicates that for Pe > 10, thermal diffusion has little effect on the sta- 
bility of KH modes whose horizontal wavelength is comparable to the vertical thickness 
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of the shear layer; the usual adiabatic Richardson criterion Ri = N'^/{dU/dz)'^ > 1/4 
is sufficient for the stability of these short-wavelength modes. For typical conditions 
on the day sides of hot Jupiters, Ri < 1/4 implies transsonic or supersonic winds, 
presuming that the shear layer is no narrower than a pressure scale height. 

4. Nevertheless, at any finite Pe and Ri, inviscid flows with an inflection point in the 
shear profile are unstable at sufficiently long horizontal wavelengths, at least when 
the flow is confined within a channel of finite vertical extent. We doubt that these 
long-wavelength KH modes are important for hot- Jupiter winds because their growth 
rates are very small (inversely proportional to wavelength) and because the planetary 
circumference is not large enough compared to Hp to accommodate them at Pe > 10^. 

5. We estimate a turbulent diffusivity for composition in hot Jupiter atmospheres. For 
Pe ~ 100, we find Kt/csHp ^ 10""^, depending on the strength of the horizontal 
forcing. However, these values pertain to regions of maximum vertical shear; turbulent 
diffusivities may be much smaller elsewhere. 

6. Nonlinear, two-dimensional (azimuth and depth), horizontally forced local simula- 
tions saturate at time-averaged Mach numbers ~ 2 when the Reynolds number Re = 
CsHp/u > 10'^. Given the widths of our shear layers (~ 2Hp) and background stratifi- 
cation, this corresponds to minimum Richardson numbers slightly below the adiabatic 
critical value of 1/4. In this regime, which is relevant to hot Jupiters, we find that the 
saturated wind speed depends only weakly on Pe and on the strength of the forcing. 

7. Dissipation in the high-i?e regime is due mainly to the work expended to overturn 
the stratification rather than shocks, though weak shocks are present. It is highly 
intermittent, at least in our two-dimensional calculations, and appears to be triggered 
by a recurrent linear instability of Kelvin-Helmholtz type. 

8. At Re < 10'^ or at low numerical resolution (significantly fewer than 10 cells per 
pressure scale height), multiple sequences of solutions appear in overlapping ranges of 
Reynolds number, differing in their wind speed and predominant dissipation mecha- 
nisms. The fastest have Mach numbers at least twice as large as what we believe to 
be the correct inviscid values for our forcing strengths. 

The last three points suggest that some global simulations of hot Jupiters may have 
overestimated wind speeds due to incomplete physics (e.g. neglect of vertical accelerations 
and sound waves) and/or insufficient spatial resolution. Pending sufficient computational 
resources to resolve the physical dissipation mechanisms directly in global simulations, we 
suggest that artificial dissipative terms be added to the equations of motion so as to prevent 
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strongly supersonic shears and small Richardson numbers. The appropriate form for these 
terms will depend upon the equation sets and algorithms used, but presumably it is desirable 
that they conserve energy and momentum, and probably that they have a sharp threshold 
in Mach number or Richardson number. 

We thank Jim Stone for help with ZEUS. This work was supported in part by the 
National Science foundation under grant AST-0707373. 
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A. Exact linear solutions 



We provide here four limits of eqns I14al and I14bl for which there analytic solutions. We 
test our code against these solutions. 

Taking A^^ = 0, i^ — ?■ 0, eqns I14a1 and I14bl admit solutions satisfying the inviscid pertur- 
bation equation with vanishing temperature perturbation, 



'\z 



KU" 



a 



+ k^ vu, 



Oi = 0, 



(Al) 
(A2) 



in which k = {k^ + kyY^^. For a hyperbolic tangent velocity profile U = tanh(2;), the inviscid 
perturbation equation for the marginal u = mode becomes 



'Iz 



- (2sech^2; — k"^) vi^. 
A solution to this equation that vanishes as 2; — t- ±00 is 

Viz = sechz , k'^ = 1. 



(A3) 



(A4) 



In the limit 



eqns I14al and Il4bl become 



[/ = 0, z/ -> 0, iV^ = 0, 



v'L = ik')viz - -k'Oi, 



e'l 



%UJ 



X 



+ k^ 01, 



which for Re{uj) = and Im{u) < —{k'^)x admit the decaying mode solutions 

61 = A cos 



■ \ 1/2 
-k^ + — ] z + 

xj 



Viz 



.Xk'^ 
A — — cos 



u^ 



-k^ + — 

X 



1/2 



z + 



Be^' + Ce^ 



(A5) 

(A6) 

(AT) 



(A^ 



(A9) 



In the limit x ~^ 00 with conditions IA5t we obtain solutions to eqns IA6I and IA7I of the 



form 



kz 



9i = Ae"' + Be 



~kz 



(AlO) 



-33- 



v^, = — [Ae^' + Be-^' + 2kz {Ae^' - Be"'"')) 



+06"' + De-^'. 



kz 



(Aii: 



Finally, in the limit 



eqns I14al and I14bl combine to give 



[/ = 0,x^ 0,1/^0, 



For uj real and N'^/u'^ > 1, eqn lA13l has solutions of the form 



Viz = A COS 



N^ 



k{ — -1 1 z + 



1/2 



&i = A cos 



N^ 



A; ( — - 1 1 z + 



1/2 



(A12) 



(A13) 



(AM) 
(A15) 



B. Long- wavelength modes 

Here we demonstrate analytically for the Boussinesq equations that there exist long- 
wavelength instabilities at arbitrarily large Richardson number if the Peclet number is finite. 
It is helpful to adopt dimensionless variables. Supposing U{z) ^ ±Uo at \z\ ^ H, where H 
is the width of the shear later, we define these variables as follows: 



2* = z/H 
w*(z*) = viz{z)/Uo 



U4z,) = U{z)/Uo 
9i 



p* 



mH 



K = Hjkl + kl 



Re = UqH/u 



cos a = kx/ \ k1 + k. 
Pe = UoH/x 



c* = u]/{kJJo) 
Ri = (NH/Uo)'- 



Equations (I14al) & fll4b|) become, with primes denoting d/dz 



U" 



(c* — f/*) cos a 



i^-kl 



v^ 



JXl K:kt7^ 



k^.Re \dzl 
O'l + [iPe K cos a{c^-U^) - k'f\ 9^ = Pev^. 



(Bl) 

(B2a) 
(B2b) 



We take the inviscid, long wavelength limit in such a way that k^^Re — )■ oo while A;* — ?■ 0^, 
with Pe and Ri regarded as finite and nonzero. The wavelength 27r|fc|^^ is not infinite, but 
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rather ^ 27ri7, so that we may neglect the exphcit appearances of /c* in the equations above 
yet regard c* as finite. If c^ should turn out to have an imaginary part, then the growth rate 
Im(a;) — > A;^[/oIm(c,,) as k^ — )■ 0. The potential temperature Q^ decouples: 

<' + — ^t^* = 0, (B3a) 

e'l = Pev,. (B3b) 

The first of equations flB3P is the long-wavelength limit of the inviscid, unstratified Kelvin- 
Helmholtz problem with a continuous shear profile t/=f(z*). Given boundary conditions 
w*(±z*,niax) = for some z*,max that is finite and sufficiently large, and presuming that 
f/* is an odd function of z^ so that f/"(0) = 0, then there is a solution to ( 1B3|) for which c* 
has a positive imaginary part. In fact, although it doesn't satisfy the boundary conditions, 
!>* oc [/* — c* is a solution to eq. (lB3aP for any c*, and a second solution v^ = f{z^){U^ — c*) 
can be found by solving a first-order equation for f\z^), followed by a quadrature. Requiring 
that this solution vanish at both ±2*^max determines c* via 

"^^^ 0. (B4) 



•2*, max 



[f/*(-2*) - C*]2 



Since [/* is odd and ~ ±1 at large |z*|, it follows easily that c* ^ ±z when 2;*, max ^ 1- 

This demonstrates the existence of inviscid instability when |A;|^^ ^ -2* max ^1- If 
2^*,max ^ ^* ! however, then a separate analysis would be needed, because the solution of 
eq. ( lB3bl) might yield 6** ~ Pezl^^v*, and then the righthand side of eq. ( lB2al) could be 
0(1) or larger. We have not done this analysis because the numerical evidence indicates that 
the long-wavelength growth rates are too slow to be important for our applications. 



